Towards an assessment of the accuracy of density functional theory for first principles 

simulations of water II 



O 

o 

(N 



(N 



o 

CO 



C 

o 

(N 
> 

in 
in 
o 
^i- 
o 

-I— > 



i 

c 

o 
o 



X 
S3 



Eric Schwegler, Jeffrey C. Grossman, Frangois Gygi, and Giulia Galli 
Lawrence Livermore National Laboratory, P.O. Box 808, Livermore, California 84559 

(Dated: February 6, 2008) 

A series of 20 ps ab initio molecular dynamics simulations of water at ambient density and 
temperatures ranging from 300 to 450 K are presented. Car-Parrinello (CP) and Born-Oppenheimer 
(BO) molecular dynamics techniques are compared for systems containing 54 and 64 water molecules. 
At 300 K, excellent agreement is found between radial distribution functions (RDFs) obtained with 
BO and CP dynamics, provided an appropriately small value of the fictitious mass parameter is used 
in the CP simulation. However, we find that the diffusion coefficients computed from CP dynamics 
are approximately two times larger than those obtained with BO simulations for T>400 K, where 
statistically meaningful comparisons can be made. Overall, both BO and CP dynamics at 300 K 
yield overstructured RDFs and slow diffusion as compared to experiment. In order to understand 
these discrepancies, the effect of proton quantum motion is investigated with the use of empirical 
interaction potentials. We find that proton quantum effects may have a larger impact than previously 
thought on structure and diffusion of the liquid. 



I. INTRODUCTION 

Recently, a number of studies have focused on un- 
derstanding in detail the level of accuracy that can be 
achieved in ab initio molecular dynamics simulations of 
wateni* 2 * 3 -. In particular, it was found in both Rcf. 2 and 3 
that in the absence of proton quantum effects, density 
functional theory (DFT) within widely used generalized 
gradient approximations (PBE 4 and BLYP^&) leads to 
significantly overstructured radial distribution functions 
(RDFs) and slower diffusion than experiment, at ambient 
temperature and pressure. 

In our previous work 3 -, within the Car-Parrinello 
(CP) method^, we found that using either the Perdew- 
Burke-Ernzerhof (PBE) 4 or the Becke-Lee- Yang-Parr 
(BLYP)^ functional leads to very small differences in 
the structural properties of the room temperature liquid. 
Size effects, although not fully negligible when using 32 
molecule cells, were found to be rather small. We also 
found that there is a wide range of ratios {n/M) between 
the fictitious electronic mass (/i) and the smallest ionic 
mass of the system (either M ~ 1836 au for hydrogen 
or M ~ 3672 au for deuterium atoms) for which the 
electronic ground state is accurately described in micro- 
canonical CP simulations (n/M < 1/5). However, care 
must be exercised not to carry out simulations outside 
this range, where structural properties may artificially 
depend on /i. Specifically, our results showed that the use 
of fi/M ~ 1/3 leads to an artificial softening of the pair 
correlation functions, which are in fortuitous agreement 
with experiment. In the case of an accurate description 
of the electronic ground state, and in the absence of pro- 
ton quantum effects, the oxygen-oxygen RDF is found to 
be over-structured compared with that obtained in neu- 
tron scattering 8 and x-ray diffraction experiments^ and 
the computed diffusion coefficients are 10 to 20 times 
smaller than experiment. In our previous work 3 -, we also 
examined the amount of sampling needed to obtain sta- 
tistically independent data points for calculated proper- 



ties, and showed that for segments of the oxygen-oxygen 
RDF the data correlation time is quite long; time inter- 
vals as large as ~10-20 ps are required to obtain a single 
statistically independent data point (correlation time is 
not to be confused with equilibration time, as we will 
discuss later in the text.) 

Although a number of important technical parame- 
ters were tested in Ref£, several questions remain unan- 
swered. First, it is important to fully test the convergence 
of the CP algorithm as the value of the fictitious mass ap- 
proaches zero and to investigate whether possible sources 
of inaccuracies, such as those pointed out by Tangney and 
Scandoloi^ persist even for small fi/M ratios (/J./M < 
1/5). In order to address this issue, we present here a se- 
ries of ab initio MD simulations at different temperatures, 
using both the CP method and minimizing the Kohn- 
Sham energy functional at each MD time step, while 
keeping all the other parameters of the two calculations 
fixed. We refer to this direct minimization technique as 
Born-Oppenheimer (BO) MD. Our results show an excel- 
lent agreement between RDFs obtained in CP and BO 
simulations, provided an appropriately small value of the 
fictitious electronic mass is used (ft/M < 1/5), in agree- 
ment with previous BO simulations 2 . However, diffusion 
coefficients obtained within BO and CP simulations ap- 
pear to differ by about a factor of two at temperatures 
(T) of ~ 400K. At this temperature, a statistically mean- 
ingful comparison is possible using the simulation times 
employed here (about 20 ps). Unfortunately, at lower 
temperatures such a quantitative comparison will require 
longer simulation times. 

In addition to comparing BO and CP simulations, we 
have carried out a detailed analysis of the effect of proton 
quantum motion on computed structural and dynami- 
cal properties, by using the results of extended classical 
simulations with empirical inter-atomic forces. Our find- 
ings indicate that the inclusion of proton quantum ef- 
fects would bring the DFT results obtained here in much 
better agreement with experiment, for both RDFs and 
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diffusion coefficients, and that their inclusion is crucial 
to accurately describe the properties of water at ambient 
conditions. In agreement with most previous studies, we 
find that proton quantum effects tend to give softer RDFs 
and faster diffusion than classical (Newtonian) simula- 
tions. We also find that the effect of proton quantum 
motion on calculated RDFs and diffusion coefficients at 
300 K is the same as that obtained by increasing the ionic 
temperature by about 50 or 100 K at constant density, 
when using a rigid or flexible water model, respectively. 

The rest of the paper is organized as follows. ScctionlTTl 
contains a description of the methods used in both ab ini- 
tio and classical MD simulations of water and Section ITTT1 
describes our results in detail. In particular, in Section 
IIII Al wc report a comparison between CP and BO simu- 
lations, and in Section llll Bl we discuss possible sources of 
the apparent disagreement with experiment, namely the 
absence of proton quantum effects and the GGA func- 
tional (PBE) adopted here. The effect of temperature 
variations between 300 and 450 K on water at ambient 
density is analyzed in Section fill CI Finally, Section llVl 
contains our conclusions. 



II. METHODS 

We have performed a series of molecular dynamics 
simulationsii of liquid water with both ab initio and 
classical potentials. The ab initio molecular dynamics 
simulations are based on density functional theory with 
the PBE4 generalized gradient approximation. Electron- 
ion interactions are treated with norm conserving non- 
local pseudopotentials of the Hamann typ o 12 i 13 , and with 
the Kohn-Sham orbitals and charge density expanded in 
plane waves up to a kinetic energy cutoff of 85 and 340 
Ry, respectively. In each simulation, the starting config- 
uration was generated by performing 500 ps of classical 
molecular dynamics with the TIP5P potential 14 , followed 
by 3 ps of ab initio molecular dynamics with a weakly 
coupled velocity scaling thermostat set at the given tar- 
get temperature. The thermostat was then removed and 
all reported statistics were collected under constant en- 
ergy conditions. 

For the Born-Oppcnhcimer (BO) simulations, the sys- 
tem consists of 64 water molecules in a cubic box of 
length 12.43 A with periodic boundary conditions. At 
each molecular dynamics timestep, the wavefunction was 
relaxed to the ground state by preconditioned steepest 
descent with Anderson acceleration 1 ^. For a timestep of 
0.24 fs, 12 electronic iterations were found to be suffi- 
cient to converge the total energy to within 10~ 9 a.u. at 
each molecular dynamics iteration. Under these simula- 
tions conditions, the drift in the total energy was found 
to be equivalent to a 0.27 K/ps increase in the system 
temperature. 

In the Car-Parrinello (CP) simulations, a slightly 
smaller system size of 54 water molecules in a cubic box 
of length 11.7416 A was used. The electronic degrees of 



TABLE I: Geometry and potential parameters defined as in 
Refspi 4 *^ for the TIP5P water models. 





TIP5P 


TIP5P(PIMC) 


Qir(e) 


0.241 


0.251 


q L {e) 


-0.241 


-0.251 


°o (A) 


3.12 


3.12 


to (kcal/mol) 


0.16 


0.16 


r H (A) 


0.9572 


0.9572 


roh (A) 


0.70 


0.70 


6hoh (deg) 


104.52 


104.52 


Olol (deg) 


109.47 


109.47 



freedom were propagated with a fictitious mass of ^t=340 
a.u., and with a molecular dynamics timestep of 0.07 fs. 

All classical simulations were performed with the Tin- 
ker simulation packaged under constant volume and 
temperature conditions for systems of 64 rigid water 
molecules. A molecular dynamics timestep of 0.5 fs was 
used and the simulation temperature was maintained 
with a weakly coupled Berendsen-type thermostat. 
Long-range electrostatic interactions were included with 
the particle-mesh Ewald method 19 . Both the TIPSRii 
and the TIP5P(PIMC)i& potentials, which are defined in 
Tabled have been examined. 



III. RESULTS 
A. Comparison of CP and BO dynamics at 300 K 

In Ref£, we have shown that in CP simulations of wa- 
ter at ambient conditions, one can vary the fictitious mass 
parameter, (i, within a suitable range without signifi- 
cantly altering the results of the simulations. In particu- 
lar, it was found that when [i/M < 1/5, where M is the 
mass of the H or D atom in the simulation, the average 
structural properties and computed diffusion coefficients 
appear to be insensitive to the choice of fi. This is in 
contrast to simulations performed with large values of /j, 
that lead to significant changes in both the structure and 
diffusion of the liquid. As discussed in Ref£, the reason 
for the dramatic changes observed as n/M is increased 
from 1/5 are due to a direct overlap of the fictitious elec- 
tronic spectra with the highest frequency ionic spectra. 
In other words, as long as the fictitious electronic degrees 
of freedom can be adiabatically decoupled from the ionic 
degrees of freedom, the Car-Parrinello technique yields 
accurate results. However, as argued by Tangney et alMl, 
exact adiabatic separation in a Car-Parrinello simulation 
is never actually achieved, even for relatively small val- 
ues of /i. The reason for this is because in addition to 
the high frequency components in the ionic spectra that 
can lead to direct coupling to fictitious electronic degrees 
of freedom, there are additional /i-dependent errors as- 
sociated with the low-frequency components of the ionic 
motion. The impact of a particular choice of /i can be 
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FIG. 1: Comparison of the oxygen-oxygen radial distribution 
functions obtained in BO and CP simulations at a temper- 
ature of 300 K. The solid line corresponds to CP dynamics, 
the dashed line to BO dynamics, and the dot-dashed line to 
experiment. 



FIG. 2: Comparison of the oxygen-hydrogen radial distribu- 
tion functions obtained in BO and CP simulations at a tem- 
perature of 300 K. The solid line corresponds to CP dynamics, 
the dashed line to BO dynamics, and the dot-dashed line to 
experiment. 



assessed in a direct fashion by comparing Car-Parrinello 
simulations to so-called Born-Oppcnhcimcr (BO) simula- 
tions of water, where the wavefunctions are converged to 
the ground state at each molecular dynamics time step. 

In Figs. to the oxygen-oxygen, oxygen-hydrogen 
and hydrogen-hydrogen radial distribution functions 
(RDFs)2fi obtained with CP (Simulation 1) and BO (Sim- 
ulation 4) ab initio MD at a temperature of approxi- 
mately 300 K are shown along with the corresponding 
distribution functions obtained by experiment^. As can 
be seen, both the BO and CP simulations of water at 
300 K lead to an overstructured liquid as compared to 
experiment. In addition, these results confirm the previ- 
ous observation that, as long as \x is chosen appropriately, 
CP simulations of water at ambient conditions can give 
a consistent description of the liquid RDFsi 

The large amount of overstructure seen in Figs. ^ to 
13 appears in the simulation within a relatively short 
timescale (~2 ps), and can be readily observed, even 
when starting the simulation from a classical theoret- 
ical model or a higher simulation temperature, which 
may yield a less structured g(r). This relatively short 
timescale should not be confused with the much longer 
timescale required to obtain statistically independent 
data points for a given RDF 3 , which will be discussed 
further in Section Till B 31 

In addition to the average structural properties of the 
liquid, we have computed the diffusion coefficient using 
the Einstein relation: 

6£ = t lim K^W-r^O)! 2 }. (1) 

To evaluate Eq.Q] the mean square displacement (MSD) 
of the oxygen atoms as a function of time was determined 
by averaging over the trajectories with multiple starting 
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FIG. 3: Comparison of the hydrogen-hydrogen radial distri- 
bution functions obtained in BO and CP simulations at a 
temperature of 300 K. The solid line corresponds to CP dy- 
namics, the dashed line to BO dynamics, and the dot-dashed 
line to experiment-. 



configurations that are evenly spaced by 4 fs. The slope 
of the resulting MSD was determined in the range of 1 
to 10 ps and used to compute the diffusion coefficient, 
D. As shown in Table ITU the computed diffusion coeffi- 
cients obtained in this manner for both the CP and BO 
simulations at 300 K are much smaller (by at least a fac- 
tor of ten) than D measured experimentally at the same 
temperature^!. Although the difference between the sim- 
ulated and measured diffusion coefficients are statistically 
meaningful, the difference between the CP and BO dif- 
fusion coefficients at 300 K are not. Indeed, to make 
a quantitative comparison between diffusion coefficients 
that are in the range 10" 7 to 10" 6 cm / sec, orders of 
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TABLE II: For each simulation performed in this work, the type of simulation (Car-Parrinello - CP or Born-Oppenheimer - 
BO, number of molecules N, total simulation time (ps), average temperature (K), diffusion coefficient D (cm 2 /s), position (A) 
and value of first maximum and minimum in the goo(r)~, average temperature (K), and average coordination (CN) of the 
water molecules are listed. All simulations are for H2O with the PBE density functional. CP simulations were carried out with 
a fictitious electron mass of 340 (1100) a.u. for flexible (rigid) water. The last row contains measured diffusion coefficients 21 
and structural data 8 for water at ambient conditions. 



Sim. 


TvDe 




N 


Time 


Tavg 


D 




P [srl r I™,....! 

H&V L Jmax\ 


max 


P [{rfr)™..'.™! 

H&V L Jrmn\ 


Tnin 


CN 


1 


CP 




54 


19.8 


296 


2.4x10" 


■e 


2.69 


3.65 


3.32 


0.37 


4.2 


2 


CP 




54 


18.8 


345 


5.0x10" 


6 


2.72 


3.21 


3.35 


0.42 


4.3 


3 


CP 




54 


22.0 


399 


2.2x10" 


■5 


2.75 


2.60 


3.41 


0.73 


4.6 


4 


BO 




64 


20.5 


306 


7.9x10" 


■7 


2.72 


3.83 


3.25 


0.33 


4.1 


5 


BO 




64 


20.7 


349 


3.3x10" 


6 


2.72 


3.49 


3.30 


0.40 


4.2 


6 


BO 




64 


18.6 


393 


1.2x10" 


■5 


2.73 


3.10 


3.40 


0.56 


4.6 


7 


BO 




64 


10.5 


442 


3.9x10" 


■5 


2.75 


2.63 


3.44 


0.74 


4.8 


8 


CP-Rigid 




54 


24.5 


315 


1.3x10" 


■5 


2.75 


2.92 


3.41 


0.61 


4.6 


9 


CP-Rigid 




54 


26.3 


345 


3.3x10" 


■5 


2.75 


2.61 


3.41 


0.77 


4.7 


10 


TIP5P(PIMC),P= 


:1 


64 


500.0 


300 


2.6x10" 


6 


2.69 


3.61 


3.29 


0.43 


4.1 


11 


TIP5P(PIMC),P= 


■1 


64 


500.0 


350 


2.8x10" 


■5 


2.72 


2.84 


3.38 


0.78 


4.6 


12 


TIP5P(PIMC),P= 


=5 


216 




300 






2.73 


2.76 


3.44 


0.77 






Experiment 








298 


2.3x10" 


■5 


2.73 


2.75 


3.36 


0.78 


4.7 



magnitude longer simulation times than those used here would be required. 



B. Comparison of ab-initio MD to experiment and 
analysis of discrepancies 

There are a number of possible explanations for 
the observed overstructure and slow diffusion of our 
DFT/GGA-based molecular dynamics simulations of wa- 
ter. In the following, we concentrate on three of the most 
likely causes. The first is related to possible inaccuracies 
of the PBE functional in describing hydrogen bonding. 
The second is related to the neglect of proton quantum 
effects. The third is related to the timescales used in the 
ab initio simulations. 



1. Accuracy of DFT/GGA for the water monomer and 
dimer 



TABLE III: Selected properties of the water monomer and 
dimer. All reported distances are in A. DM refers to the 
dipole moment in Debye, cti ao is the isotropic polarizability 
in A 3 , OHO is the intermolecular hydrogen bond angle, and 
D e is the binding energy of the water dimer in kcal/mol. 







Monomer 






Dimer 




Method 


roH 


IHOH 


DM 




roo 


/OHO 


D e 


LDA 


0.970 


104.9 


1.87 


1.52 


2.71 


172.8 


9.02 


PBE 


0.970 


104.2 


1.81 


1.55 


2.90 


173.7 


5.10 


BLYP 


0.972 


104.5 


1.80 


1.54 


2.95 


171.6 


4.18 


PBE1PBE 


0.959 


104.9 


1.85 


1.42 


2.90 


174.1 


4.90 


B3LYP 


0.962 


105.1 


1.85 


1.45 


2.92 


172.1 


4.57 


CCSD(T) 


0.959° 


104.2 a 






2.91 6 


174.5 6 


5.02 6 


Exp. 


0.957 c 


104.5 C 


1.85 d 


1.43 e 


2.98 / 


174.0 / 


5.44 s 



6 Ref^i. c Ref^. d Ref^&. e Ref^£. ^Ref^. 9Ref-' 
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ill order to better quantify the accuracy of the 
DFT/GGA functional used here, it is useful to first exam- 
ine the computed properties of the water monomer and 
dimer (for the interested reader, a more comprehensive 
review of the water monomer and dimer properties as 
computed within DFT was recently published in Ref. 22 ). 

As shown in Table IIIII we have examined the use of 
a number of DFT exchange-correlation functionals, in- 
cluding LDA22i2i, GGA (PBBl and BLYR££), and hy- 
brid (PBEIPBB^ and B3LYRS) functionals. Also in- 
cluded in Table HTU are the results of several CCSD(T) 
calculationsSiSi as well as the relevant experimental 
measurementa 25 ! 26 ! 27 ! 28 ! 29 . The CCSD(T) results repre- 
sent the most accurate quantum chemistry computations 
of the water monomer and dimer to date. For all of the 



DFT-based calculations reported in Table IIIII we have 
used the electronic structure program Gaussian 03^4 with 
the aug-cc-pVTZ basis sel^. 

For the computed bond distances and angles, the ma- 
jority of the tested functionals perform reasonably well. 
The most notable exception is the oxygen-oxygen dis- 
tance in the LDA water dimer calculation, which is much 
too short. Although the computed dipole moments of the 
water monomer are also in good agreement with each 
other, it is interesting to note that the isotropic polar- 
izabilities are systematically too large at the LDA and 
GGA level of theory. This is due to the well-known 
tendency of local DFT functionals to underestimate the 
HOMO-LUMO gap of molecular systems, which in turn 
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leads to an overestimation of the molecular polarizability. 
As more accurate, hybrid functionals are considered, this 
general tendency seems to disappear. In Table llTTl we also 
show computed binding energies. As expected, the LDA 
binding energy is much too large. The remainder of the 
DFT functional binding energies are in relatively good 
agreement with each other, as well as with the CCSD(T) 
calculation and the experimental measurement. 

Given the differences in the water dimer binding en- 
ergies found with BLYP and PBE (the BLYP binding 
energy is ~1 kcal/mol smaller than PBE), one might 
think that the PBE functional would give more struc- 
tured RDFs and slower diffusion of the liquid, as com- 
pared to BLYP. In our previous work we could not resolve 
any statistically significant difference between PBE and 
BLYP results. However, comparison of the present CP 
results with those of Chen et al^ - which were obtained 
with a sufficiently small fictitious mass (/j,=400 a.u.) and 
the BLYP functional - seem to show differences which 
would be consistent with the idea that BLYP yields less 
structured RDFs than PBE. At this point, it is difficult to 
draw any firm conclusions, as the time scale of both the 
present work and the simulations of Ref.— are probably 
too small to statistically resolve differences. 

We note that there is a very good agreement between 
the PBE and the CCSD(T) binding energies, and overall 
the water monomer and dimer properties listed in Table 
IIIII appear to be well reproduced by the PBE functional. 
Although these zero temperature, gas phase results do 
not necessarily insure an equally good performance of 
the functional in the liquid state, they are very promis- 
ing. We also note that the PBE functional appears to 
be quite accurate in reproducing the sublimation energy, 
lattice constant and bulk modulus of ice-Ih-I. It would be 
interesting to explore the accuracy of other functionals, 
along the lines proposed in Rcf. 2,3 — , especially of hybrid 
functionals, but this falls beyond the scope of the present 
investigation. 

2. Proton quantum effects 

In addition to possible inaccuracies in the DFT/PBE 
treatment of water, it is plausible that some amount of 
the observed overstructure is due to the neglect of the 
quantum motion of the hydrogen atoms. As pointed 
out in Ref. 39 , at 300 K, k B T ~ 200 cm" 1 , whereas the 
high frequency intramolecular modes in water range from 
^1000 to 3500 cm -1 . Therefore, in the real quantum 
system, the amount of thermal energy available to ex- 
cite vibrational modes is much smaller than the lowest 
possible intramolecular vibrational excitations, 

Tilo > k B T. (2) 

In other words, the quantum system will be restricted to 
its vibrational ground state at a temperature of 300 K, 
and the quantum and the classical system will be qual- 
itatively different in terms of the distribution of vibra- 
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FIG. 4: The oxygen-oxygen radial distribution function ob- 
tained in TIP5P(PIMC) simulations (see text) of water at 
a density of 0.997 g/cc, and a temperature of 300 K. The 
solid line corresponds to a simulation performed with P=l 
(no path integral sampling), the dashed line to P=5 (taken 
from Ref^-), and the dotted line to experiment-. 



tional energy. When using empirical potentials with in- 
tramolecular flexibility, quantum effects can, to some ex- 
tent, be implicitly included in the potential parameter- 
ization. However, in the case of ab initio based meth- 
ods, the possibility of implicitly accounting for quantum 
effects by renormalizing the interaction potential is, by 
definition, not an option. 

In principle, the effect of the proton quantum motion 
can be accounted for with path-integral (PI) methods^. 
To date, combining PI sampling of ionic trajectories with 
empirical interaction potentials has been found to lead 
to an overall softening of the RDFs of water at ambient 
conditionsiS^i^Siiiii. However, the majority of these 
simulations have been performed with empirical poten- 
tials that already include some amount of quantum ef- 
fects. In particular, it is common to use classical molec- 
ular dynamics or Monte Carlo sampling in the parame- 
terization of an interatomic potential, with the goal of 
fitting experimental data as well as possible. Within 
this procedure, agreement with experiment is achieved 
by renormalizing interatomic interactions so as to in- 
clude proton quantum effects, which are not explicitly 
treated by classical molecular dynamics or Monte Carlo. 
There are, however, important exceptions to this general 
procedure. In particular, in Ref. 1 —, the TIP5P potential 
was reparameterized with PI Monte Carlo sampling in- 
stead of classical Monte Carlo sampling. The resulting 
potential, which we will refer to here as TIP5P(PIMC), 
is nearly identical to the original TIP5P potential ex- 
cept for the magnitude of the charges on the hydrogen 
and lone-pair sites, which have been slightly increased 
(see Tabled}. The TIP5P(PIMC) potential opens up the 
interesting possibility of directly examining, in the con- 
text of a rigid water model, the extent to which quan- 
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FIG. 5: The probability distribution of water molecule dipole 
moments obtained in CP simulations of water at a temper- 
ature of 300 K. The solid line corresponds to a simulation 
performed with the fictitious mass /i=340 au, and the dashed 
line to a simulation with /i=760 au. 



sical ion dynamics to an average value of 3.19 Debye 
when PI sampling is included. Although we have used 
a different GGA functional (PBE instead of the BLYP 
functional used in Ref. 45 ), it is interesting to note that 
in our CP simulations we find a dipole moment of 3.24 
Debye when a fictitious mass of /i=320 au is used (a 
value of /z=500 au was reported in Ref. 45 ). As shown 
in Fig. [5] when a larger value of the ficticious mass is 
used (/x=740 au) the average dipole moment decreases to 
3.06 Debye. An average value of 3.24 Debye found here 
is not consistent with the results reported in Ref^ 5 -, and 
opens up the possibility that the role of quantum effects 
in path integral/DFT-based simulations of water may not 
be completely understood. Indeed, simulations using PI 
sampling with a flexible water model are known to be 
quite difficult to converge 16 . In this regard, additional 
investigations, which are beyond the scope of this work, 
are most likely needed where the convergence of the var- 
ious technical parameters involved in PI simulations are 
examined in detail. 



turn effects can soften the RDF of water. As shown 
in Fig. E| when the TIP5P(PIMC) potential is used in 
a classical molecular dynamics simulation with P=l (P 
refers to the number of beads used to represent the path 
integrals), the resulting oxygen-oxygen RDF becomes 
severely overstructured as compared to a PI Monte Carlo 
calculation with the TIP5P(PIMC) potential where P=5 
is used (taken from Refii 6 -). The amount of overstruc- 
ture seen in Fig. 0] is nearly identical to the BO and CP 
RDFs presented earlier in Fig. ^ I n addition, in the 
TIP5P(PIMC), P=l simulation the diffusion coefficient 
of the liquid is 2.2 x 10 _6 cm 2 /s, which is approximately 
ten times smaller than the experimental value; this value 
is also in qualitative agreement with the corresponding 
BO and CP simulations performed at 300 K (see Table 
ITT}) . These results indicate that including proton quan- 
tum effects may bring the ab initio results presented here 
into significantly better agreement with experiment and 
that neglect of these effects is probably responsible for 
the majority of the discrepancy between experiment and 
DFT results. 

Along with the existing PI studies based on empiri- 
cal water potentialsi 6 ^^^^^, Chen, et alJ& recently 
demonstrated the feasibility of directly performing PI 
sampling and CP simulations of water. However, instead 
of the expected softening of the liquid structure, Chen, 
et found the rather surprising result that quantum 
effects may actually increase the structure of water. The 
reason for this increase was attributed to an enhancement 
of the individual water molecule's dipole moments when 
proton quantum effects are included in the calculations. 
In particular, an analysis based on maximally localized 
Wannier functions^ 6 " was used to demonstrate that the 
average dipole moment of water molecules increases from 
3.08 Debye in their CP simulation performed with clas- 



3. Timescales 

As discussed in Ref. 3 , another approximation that 
needs to be considered is the length of the time inter- 
val over which averages, such as the RDFs, are taken. 
In the case of ab initio molecular dynamics simulations, 
which are usually limited to relatively short timescales 
(10 to 20 ps), this can be a particularly significant effect. 
In this regard, it is important to determine the time in- 
terval outside which neighboring simulation data points 
are no longer correlated, and thus statistically significant. 
In Ref. 3 a method based on a combination of re-blocking 
and autocorrelation function techniques was used for de- 
termining this time interval along with an estimate of 
the variance for individual bins of a RDF. Because this 
method requires simulation trajectories of at least 1 to 
2 orders of magnitude longer than the actual correla- 
tion time, a 2 ns classical MD simulation of 32 TIP5P 
water molecules under ambient conditions was used for 
the analysis. Based on this trajectory, correlation times 
as long as 10 ps were observed for the oxygen-oxygen 
RDF. However, we note that the diffusion coefficient for 
32 TIP5P water molecules under ambient conditions is 
~ 2.1 x 10 _5 cm 2 /s, which is at least 10 times larger 
than what is found in our ab initio BO and CP sim- 
ulations at 300 K. Therefore, in order to have an esti- 
mate of the RDF correlation time and variance relevant 
to our 300 K BO and CP simulations, it is more appro- 
priate to use the TIP5P(PIMC), P=l potential for this 
analysis. In Fig. the computed correlation time for 
the oxygen-oxygen RDF bins is shown for both TIP5P 
and the TIP5P(PIMC), P=l simulations of ambient wa- 
ter. As can be seen near the regions of the first peak in 
30o(^)j the correlation times increase from ~ 10 ps with 
the TIP5P potential to - 20 ps for the TIP5P(PIMC) 
potential. Although there is a considerable amount of 
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FIG. 6: The oxygen-oxygen RDF correlation times com- 
puted from classical molecular dynamics simulations using 
the TIP5P (solid line) and TIP5P(PIMC), P=l (dashed line) 
potentials at ambient conditions. The points correspond to 
the correlation time for each bin that was used to collect the 
distribution function. 




FIG. 7: Oxygen-oxygen RDF for the TIP5P(PIMC), P=l 
classical potential. The error bars correspond to the square 
root of the variance. 

scatter in the correlation times shown in Fig.^J this anal- 
ysis does provide a useful lower bound to the expected 
correlation times of the radial distribution function bins. 

In Fig. [7| an estimate of the variance in the oxygen- 
oxygen RDF bins is shown for a block of length equal to 
the TIP5P(PIMC) correlation times displayed in Fig.© 
As discussed in Ref£, the error bars shown in Fig. [7| 
are representative of the uncertainties that should be ex- 
pected when comparing independent ~ 20 ps simulations 
of water that are characterized by a diffusion coefficient 
on the order of I0~ 6 cm 2 /s. 

In addition to uncertainties in structural properties, 
the limited timescales of a typical first-principles simula- 
tion of water can strongly influence dynamical properties 




Time (ps) 

FIG. 8: The computed diffusion coefficient as a function of 
simulation time for the TIP5P(PIMC), P=l potential at a 
temperature of 300 K and a density of 0.997 g/cc. The circles 
correspond to a simulation that was initially equilibrated with 
the TIP5P potential at 600 K and the squares to a simulation 
that was equilibrated with the TIP5P potential at 300 K (see 
text for details). 



such as computed diffusion coefficients. In order to illus- 
trate this effect we have performed two simulations of 
64 water molecules with the TIP5P(PIMC), P=l poten- 
tial that differ in how the starting configurations were 
equilibrated. In the first simulation, the initial starting 
configuration was taken from a 500 ps TIP5P simulation 
at a temperature of 300 K. The simulation was then con- 
tinued for an additional 53 ps with the TIP5P(PIMC), 
P=l potential at a temperature of 300 K. This simulation 
closely resembles the equilibration procedure that was 
used for the first-principles simulations presented here. 
In the second simulation, the initial starting configura- 
tion was taken from a 500 ps TIP5P simulation at a tem- 
perature of 600 K. The simulation was then continued, as 
before, for an additional 53 ps with the TIP5P(PIMC), 
P=l potential at a temperature of 300 K. The last 50 ps 
(the first 3 ps were discarded) of these two simulations 
were analyzed by computing the diffusion coefficient via 
Eq. ^ as a function of the simulation time at 5 ps inter- 
vals. As can be seen in Fig. [SI the specifics of the initial 
equilibration phase on the diffusion coefficients can be 
rather dramatic. In particular, an equilibration phase 
that involves annealing at a higher temperature can lead 
to a significant overestimation of the diffusion coefficient 
with respect to the infinite time limit, even for simula- 
tions as long as ~ 50ps. 



C. The effect of temperature in simulations of 
water 

As pointed out in Ref. — for the ST2 water model, and 
in Ref4£ for the SPC and TIP4P models, the main struc- 
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FIG. 9: The oxygen-oxygen radial distribution function ob- 
tained in different TIP5P(PIMC) simulations. The dotted 
line corresponds to a simulation performed at 300 K with 
P=l, the solid line to 350 K with P=l, the dashed line to 300 
K with P=5 (taken from Ref^), and the dashed-dotted line 
to experiment-. 



FIG. 10: The oxygen-oxygen radial distribution function ob- 
tained in CP simulations of rigid water. The solid line cor- 
responds to a simulation performed at 315 K, the dashed 
line to a simulation at 345 K, and the dashed-dotted line 
to experiment 8 . 



tural changes that occur when going from a classical to 
a quantum description of water at 300 K are similar to 
simply increasing the temperature of the classical sim- 
ulation by 50 K. In order to examine if similar results 
are obtained when the TIP5P(PIMC) potential is used, 
we have performed classical simulations (P=l) of wa- 
ter with the TIP5P(PIMC) potential at temperatures of 
300 and 350 K (simulations 10 and 11 in Table |n|). In 
Fig. ED the oxygen-oxygen RDF obtained in these simu- 
lations are compared to TIP5P(PIMC), P=5 at 300 Ki&, 
and to experiment at 300 K 8 . As expected from pre- 
vious investigations^^, the main structural differences 
between the classical and quantum fluid described by the 
TIP5P(PIMC) are remarkably similar to a 50 K temper- 
ature increase in the classical simulation temperature at 
constant density. In addition to the observed structural 
changes, the 50 K temperature increase also causes the 
diffusion coefficient of classical TIP5P(PIMC) water to 
increase from 2.2 x 10 _6 cm 2 / s to 3.0 x 10 _5 cm 2 /s, which 
is in much better agreement with the experimental mea- 
surement of 2.2 x 10 -5 cra 2 /s for water under ambient 
conditions^. We also note that approximately account- 
ing for quantum effects via temperature scaling is a tech- 
nique that has been used in a variety of materials other 
than wate r 47 ' 48 . 

Given the agreement of results obtained with the 
TIP5P(PIMC) potential when including quantum effects, 
and increasing the classical simulation temperature, it 
is interesting to examine if a similar decrease in struc- 
ture occurs in DFT simulations of water as a function 
of temperature. However, it is important to note that 
the TIP5P(PIMC) potential used above is based on a 
rigid water model that does not allow for intermolecular 
flexibility. Therefore, to make a direct comparison be- 



tween classical and ab initio MD, we have first examined 
the effects of increasing the temperature in CP simula- 
tions of rigid water. Our results are reported in Fig. ED 
The technical details of the simulations are the same as 
those used in simulation B of Ref A As can be seen 
by comparing Figs. and ED although the amount of 
structure in the CP rigid water simulation at a temper- 
ature of 315 K (simulation 8) is smaller than the one 
observed in the flexible water simulation (simulation 1 at 
296 K), there is still a noticeable amount of overstructure 
as compared to experiment. Also shown in Fig. 1101 are 
the results of a CP simulation of rigid water at a tem- 
perature of 345 K (simulation 9). As was seen in the 
case of the TIP5P(PIMC) potential, at a temperature of 
^350 K the oxygen-oxygen radial distribution function 
is in good agreement with the corresponding experimen- 
tal measurement at 300 K. Therefore, we expect that an 
increase in temperature has a similar effect on the liq- 
uid structure with both the TIP5P(PIMC) potential and 
DFT/rigid water. This provides indirect evidence that 
RDFs computed with TIP5P(PIMC) and DFT/rigid wa- 
ter will soften to the same extent when PI sampling is 
performed. 

In addition to CP simulations of rigid water, we have 
also investigated the effect of temperature in CP and 
BO simulations of flexible water at constant density. As 
shown in Figs. ^] to ^| for CP dynamics, and in Figs. 1141 
to llGl for BO dynamics, an increase in the simulation tem- 
perature decreases the overstructure in the radial distri- 
bution functions found in the corresponding simulations 
at 300 K. However, close inspection of Figs. \H\to ED re- 
veals that both CP and BO simulations of flexible water 
require a significantly larger temperature increase than 
the simulations of rigid water in order to bring the radial 
distribution functions into agreement with experiment at 
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FIG. 11: The oxygen-oxygen radial distribution function ob- 
tained in CP simulations of water at different temperatures. 
The solid line corresponds to a simulation performed at ~350 
K, the dashed line to ~400 K, and the dashed-dotted line to 
experimen1»S. 



FIG. 13: The hydrogen-hydrogen radial distribution function 
obtained in CP simulations of water at different temperatures. 
The solid line corresponds to a simulation performed at ~350 
K, the dashed line to ~400 K, and the dashed-dotted line to 
experiment^. 




FIG. 12: The oxygen-hydrogen radial distribution function 
obtained in CP simulations of water at different temperatures. 
The solid line corresponds to a simulation performed at ~350 
K, the dashed line to ~400 K, and the dashed-dotted line to 
experiment^. 
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FIG. 14: The oxygen-oxygen radial distribution function ob- 
tained in BO simulations of water at different temperatures. 
The solid line corresponds to a simulation performed at ~350 
K, the dashed line to ~400 K, the dotted line to ~450 K, and 
the dashed-dotted line to experiment-. 



300 K. 

In Fig.Elthe computed diffusion coefficients for differ- 
ent simulations are shown as a function of the simulation 
temperature. Also shown in Fig. 1171 are the experimen- 
tally measured diffusion coefficients of water as a function 
of temperatur e - 21 ! 50 . Consistent with the oxygen-oxygen 
radial distribution functions shown in Fig. El the tem- 
perature dependence of the CP rigid water diffusion coef- 
ficient is similar to that of TIP5P(PIMC), P=l. This in- 
dicates that at approximately 350 K the simulation diffu- 
sion coefficient obtained with CP-rigid water simulations 
is equal to the experimental value at 300 K. 

In the flexible water simulations, the temperature de- 



pendence of the diffusion coefficients are quite different. 
Over the same temperature range of 300 to 350 K, the 
computed diffusion coefficients increase at a much slower 
rate than in the rigid water simulations, and tempera- 
tures of 400 K in the CP simulation and ~ 415 K in the 
BO simulations are needed in order to obtain diffusion 
coefficients comparable to experiment at 300 K. These 
findings are in qualitative agreement with previous obser- 
vations based on empirical potentials, which found that 
the use of a rigid water approximation tends to lead to 
an overall decrease in structure of the liquid and faster 
diffusion as compared to a flexible representation with 
the same intramolecular potential. 
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FIG. 15: The oxygen-hydrogen radial distribution function 
obtained in BO simulations of water at different temperatures. 
The solid line corresponds to a simulation performed at ~350 
K, the dashed line to ~400 K, the dotted line to ~450 K, and 
the dashed-dotted line to experiment. 



FIG. 17: The diffusion coefficient of water as a function of 
temperature. The diamonds correspond to the experimental 
measurements reported in Refi—, crosses to rigid water CP, 
triangles to TIP5P(PIMC), P=l, circles to flexible water CP, 
and squares to flexible water BO. The lines are simply a guide 
to the eye. 



IV. CONCLUSION 




FIG. 16: The hydrogen-hydrogen radial distribution function 
obtained in BO simulations of water at different temperatures. 
The solid line corresponds to a simulation performed at ~350 
K, the dashed line to ~400 K, the dotted line to ~450 K, and 
the dashed-dotted line to experiment-. 



We note that for all of the temperatures considered 
here, the diffusion coefficients obtained with CP dynam- 
ics (with ^=340 a.u.) appear to be approximately two 
times larger than in the corresponding BO simulations. 
As mentioned earlier in Section IIII Al and discussed in 
more detail by Tangney et aliS , this difference between 
the computed diffusion coefficients in CP and BO simu- 
lations is most likely due to a small, but still significant, 
effect of the chosen ficticious mass value in the CP sim- 
ulations. 



We have presented a series of 20 ps ab initio MD simu- 
lations of water at ambient conditions and temperatures 
ranging from 300 to 450 K, carried out using both the 
CP and BO techniques. Radial distribution functions ob- 
tained with both approaches are in excellent agreement 
with each other, provided an appropriately small ficti- 
tious mass is used in CP simulations. Differences have 
been observed in computed diffusion coefficients, those 
obtained with BO dynamics being systematically ~ 2 
times smaller than the corresponding CP values. How- 
ever, these differences are statistically significant only 
at T ~ 400 K, given the time span of our simulations 
(20 ps). These results indicate that even for relatively 
small values of the fictitious electronic mass, some cou- 
pling between the classical motion of electronic and ionic 
degrees of freedom is unavoidable in CP simulations, as 
pointed out by Tangney et ali&. Furthermore, the inac- 
curacies introduced by this coupling appear not to have 
a significant, quantitative effect on calculated structural 
properties; rather, their effect may be seen on dynamical 
properties such as the diffusion coefficient. 

We have also presented an analysis of the effects of the 
proton quantum motion on computed structural and dy- 
namical properties of water. Our results, obtained with a 
series of classical MD and PI calculations with empirical 
force fields, provide strong indications that the neglect 
of quantum effects may account for most of the discrep- 
ancy observed here between DFT-based simulations and 
experiment. We have found that including proton quan- 
tum effects is approximately equivalent to a temperature 
increase of 50 and 100 K, when using a rigid or a flex- 
ible water model, respectively, under constant volume 
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conditions. This is true for both RDFs and diffusion co- 
efficients. We note that although quantum effects most 
certainly do have an influence on the properties of water 
at ambient conditions, the exact magnitude of this effect 
is still unclear. For example, all of the simulations shown 
here have been performed under constant volume condi- 
tions that corresponds to a density of 0. 997 'g/cc. It is 
likely that this does not correspond to the exact equilib- 
rium density of water as described by the PBE functional, 
which may in turn have an effect on the precise tempera- 
tures needed to bring the RDFs and diffusion coefficients 
into agreement with experiment. 

Finally, we note that inaccuracies in the DFT descrip- 
tion of hydrogen bonding may be attributed to the choice 
of the local energy functional (PBE in this paper); the 



use of hybrid functionals that provide a better agreement 
between measured and calculated gas phase polarizabil- 
ties may improve the description of the fluid as well. 
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